Analyse géocimique de la surface de Mercure à partir des données MESSENGER
Alexandre Michaux & Beniamino Orsini
2025-10-11
1 Présentation du projet
Ce projet s’inscrit dans le cadre du projet d’analyse de données, consacré à l’étude de la planète Mercure à partir des données orbitales de la mission MESSENGER. L’objectif global est de mieux comprendre la composition chimique de la surface mercurienne et ses implications pour la structure interne et l’évolution géologique de la planète.
1.1 La mission MESSENGER
La mission MESSENGER (MErcury Surface, Space ENvironment,
GEochemistry, and Ranging)
a été lancée par la NASA le 3 août 2004 et mise en orbite autour de
Mercure le 18 mars 2011.
Elle a fourni, jusqu’à sa fin en avril 2015, un jeu de données sans
précédent sur :
- la géochimie de surface (rapports élémentaires : Mg/Si, Al/Si,
Ca/Si, Fe/Si, S/Si)
- la morphologie et la topographie (cartes globales et
régionales)
- le champ magnétique et la magnétosphère de Mercure.
Les données utilisées dans ce projet proviennent principalement de l’instrument XRS (X-Ray Spectrometer) pour la géochimie et de l’instrument MLA (Mercury Laser Altimeter) pour la topographie fournie par la publication de Nittler et al. (2020).
1.2 Objectifs du projet
Nous cherchons à analyser les variations régionales des rapports géochimiques de surface afin de répondre à une question centrale : quelles hétérogénéités de composition révèlent l’histoire magmatique et l’évolution du manteau de Mercure ?
Les ratios élémentaires (par exemple Mg/Si ou Al/Si) sont particulièrement importants, car ils renseignent sur la nature de la source mantellique et les processus de fusion qui ont produit les laves observées à la surface.
Dans ce projet, nous nous limitons aux données orbitales de surface fournies par MESSENGER. Les expériences de fusion partielle et les comparaisons avec des modèles expérimentaux (haute pression / haute température) ne seront pas encore intégrées à cette étape : l’objectif immédiat est de constituer et d’explorer un jeu de données propre et reproductible basé uniquement sur les cartes globales.
2 Préparation des données
Avant toute analyse statistique ou cartographique, il faut transformer les fichiers bruts fournis par MESSENGER en un jeu exploitable. Le travail consiste à empiler plusieurs couches de données comme dans une lasagne :
- chaque couche correspond à une variable (ex. carte globale de Mg/Si),
- toutes les couches doivent être ramenées à la même taille et au même repère spatial (ici 720×1440 en latitude–longitude),
- l’ensemble est ensuite empilé dans un cube de données 3D : latitude × longitude × variable
Ce tableau 3D constitue notre jeu de données consolidé (result_array.rds). Il permet ensuite de traverser la lasagne couche par couche (analyse univariée d’un ratio) ou de croiser les couches entre elles (analyse bivariée, cartes comparatives, corrélations).
rm(list=ls())
# ==================== PACKAGES ====================
packages <- c("raster","bmp","png","rstudioapi","abind")
installed_packages <- packages %in% installed.packages()[,"Package"]
if (any(!installed_packages)) install.packages(packages[!installed_packages])
lapply(packages, library, character.only = TRUE)## [[1]]
## [1] "raster" "sp" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[2]]
## [1] "bmp" "raster" "sp" "stats" "graphics" "grDevices"
## [7] "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "png" "bmp" "raster" "sp" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "rstudioapi" "png" "bmp" "raster" "sp"
## [6] "stats" "graphics" "grDevices" "utils" "datasets"
## [11] "methods" "base"
##
## [[5]]
## [1] "abind" "rstudioapi" "png" "bmp" "raster"
## [6] "sp" "stats" "graphics" "grDevices" "utils"
## [11] "datasets" "methods" "base"# ==================== CHEMINS ====================
base_path <- dirname(rstudioapi::getSourceEditorContext()$path)
directory_path <- file.path(base_path, "data")
result_file <- file.path(base_path, "result_array_with_uncertainty.rds")
# ==================== UTILS ====================
from03602180 <- function(mat) {
middle <- dim(mat)[2] / 2
cbind(mat[, (middle + 1):(2 * middle)], mat[, 1:middle])
}
resize_to_target <- function(mat, target_dim = c(720,1440)) {
if (is.null(dim(mat))) stop("Matrix has no dimensions.")
if (all(dim(mat) == target_dim)) return(mat)
r <- raster::raster(mat)
r_target <- raster::raster(nrows = target_dim[1], ncols = target_dim[2])
r_resized <- raster::resample(r, r_target, method = "bilinear")
as.matrix(r_resized)
}
# ==================== FONCTION DE LECTURE DES COUCHES ====================
process_files_to_3d_array <- function(directory_path, correction_factors = NULL, target_dim = c(720, 1440)) {
# -------------------------------------------------------
# Charge BMP / PNG / DAT / RDS d'un dossier,
# applique une éventuelle correction par fichier,
# redimensionne à target_dim et empile en tableau 3D.
# (CSV exclus, tous les RDS inclus)
# -------------------------------------------------------
# Liste des fichiers, en excluant explicitement CSV et TIFF
file_list <- list.files(directory_path, full.names = TRUE)
file_list <- file_list[!grepl("\\.csv$", file_list, ignore.case = TRUE)]
file_list <- file_list[!grepl("\\.tif(f)?$", file_list, ignore.case = TRUE)]
# Ne garder que les extensions gérées
file_list <- file_list[grepl("\\.(bmp|png|dat|rds)$", file_list, ignore.case = TRUE)]
n_files <- length(file_list)
if (n_files == 0) stop("Aucun fichier pris en charge (BMP/PNG/DAT/RDS) dans : ", directory_path)
layer_names <- basename(file_list)
array_3d <- array(NA_real_, dim = c(target_dim[1], target_dim[2], n_files))
apply_correction <- function(m, cf) {
if (is.null(cf)) return(m)
m * cf
}
for (i in seq_along(file_list)) {
file_path <- file_list[i]
file_name <- layer_names[i]
ext <- tolower(tools::file_ext(file_name))
message("→ Lecture de ", file_name)
if (ext == "bmp") {
# BMP : lecture brute, orientation inchangée
m <- read.bmp(file_path)
mode(m) <- "numeric"
m <- apply_correction(m, correction_factors[[file_name]])
m <- resize_to_target(m, target_dim)
array_3d[,,i] <- m
} else if (ext == "png") {
# PNG : via raster -> matrice ; flip vertical pour redresser
m <- as.matrix(raster::raster(file_path))
mode(m) <- "numeric"
m <- apply_correction(m, correction_factors[[file_name]])
m <- resize_to_target(m, target_dim)
m <- m[nrow(m):1, ] # redressement vertical
array_3d[,,i] <- m
} else if (ext == "dat") {
# DAT : structure spéciale -> recollage + (dans ta version) inversion L/C
m <- from03602180(as.matrix(read.table(file_path, header = FALSE)))
mode(m) <- "numeric"
m <- apply_correction(m, correction_factors[[file_name]])
m <- resize_to_target(m, target_dim)
array_3d[,,i] <- m
} else if (ext == "rds") {
# RDS : on prend TOUT fichier .rds
m <- readRDS(file_path)
mode(m) <- "numeric"
m <- apply_correction(m, correction_factors[[file_name]])
m <- resize_to_target(m, target_dim)
array_3d[,,i] <- m
} else {
warning("⚠️ Type non supporté/ignoré : ", file_name)
next
}
}
attr(array_3d, "layer_names") <- layer_names
array_3d
}
# ==================== FACTEURS DE CORRECTION ====================
correction_factors <- list(
"mgsi.bmp" = 0.860023 / 255.0,
"alsi.bmp" = 0.402477 / 255.0,
"ssi.bmp" = 0.161680 / 255.0,
"fesi.bmp" = 0.117737 / 255.0,
"casi.bmp" = 0.318000 / 255.0,
"mgsierr.png" = 0.860023 / 255.0,
"alsierr.png" = 0.402477 / 255.0,
"ssierr.png" = 0.161680 / 255.0,
"fesierr.png" = 0.117737 / 255.0,
"casierr.png" = 0.318000 / 255.0
)
# ==================== CONSTRUCTION DU CUBE ====================
result_array <- process_files_to_3d_array(directory_path, correction_factors)
layer_names <- attr(result_array, "layer_names")
cat("✅ Cube initial chargé :", dim(result_array)[3], "couches\n")## ✅ Cube initial chargé : 14 couches# ==================== MASQUE COMMUN (5 RAPPORTS .bmp) ====================
idx_mgsi <- grep("^mgsi\\.bmp$", layer_names, ignore.case = TRUE)
idx_alsi <- grep("^alsi\\.bmp$", layer_names, ignore.case = TRUE)
idx_casi <- grep("^casi\\.bmp$", layer_names, ignore.case = TRUE)
idx_fesi <- grep("^fesi\\.bmp$", layer_names, ignore.case = TRUE)
idx_ssi <- grep("^ssi\\.bmp$", layer_names, ignore.case = TRUE)
if (any(lengths(list(idx_mgsi,idx_alsi,idx_casi,idx_fesi,idx_ssi)) == 0)) {
stop("❌ Introuvable : au moins une des couches mgsi/alsi/casi/fesi/ssi (.bmp).")
}
# On considère qu’une valeur < 0.0001 = pixel invalide
mask_common <- (result_array[,,idx_mgsi] > 0.0001) &
(result_array[,,idx_alsi] > 0.0001) &
(result_array[,,idx_casi] > 0.0001) &
(result_array[,,idx_fesi] > 0.0001) &
(result_array[,,idx_ssi] > 0.0001)
# ==================== CUBE MASQUÉ (NA sur pixels invalides) ====================
masked_cube <- array(NA_real_, dim = dim(result_array))
for (i in 1:dim(result_array)[3]) {
L <- result_array[,,i]
L[!mask_common] <- NA_real_
masked_cube[,,i] <- L
}
# ==================== CONCATÉNER : ORIGINAL + MASQUÉ ====================
result_array_full <- abind::abind(result_array, masked_cube, along = 3)
attr(result_array_full, "layer_names") <- c(layer_names, paste0(layer_names, "_masked"))
cat("✅ Cube final :", dim(result_array_full)[3], "couches (doublé)\n")## ✅ Cube final : 28 couches (doublé)# ==================== STATISTIQUES DES COUCHES ====================
na_counts <- sapply(1:dim(result_array_full)[3], function(i) sum(is.na(result_array_full[,,i])))
layer_info <- data.frame(
Index = seq_along(attr(result_array_full, "layer_names")),
Layer = attr(result_array_full, "layer_names"),
NA_Count = na_counts
)
cat("\n=== STATISTIQUES DES COUCHES ===\n")##
## === STATISTIQUES DES COUCHES ===print(layer_info)## Index Layer NA_Count
## 1 1 alsi.bmp 0
## 2 2 alsierr.png 0
## 3 3 casi.bmp 0
## 4 4 casierr.png 0
## 5 5 DEM.RDS 0
## 6 6 DensityGrid.dat 0
## 7 7 fesi.bmp 0
## 8 8 fesierr.png 0
## 9 9 MeltGrid.dat 0
## 10 10 mgsi.bmp 0
## 11 11 mgsierr.png 0
## 12 12 ssi.bmp 0
## 13 13 ssierr.png 0
## 14 14 subregions_Nittler_Vflip.bmp 0
## 15 15 alsi.bmp_masked 369193
## 16 16 alsierr.png_masked 369193
## 17 17 casi.bmp_masked 369193
## 18 18 casierr.png_masked 369193
## 19 19 DEM.RDS_masked 369193
## 20 20 DensityGrid.dat_masked 369193
## 21 21 fesi.bmp_masked 369193
## 22 22 fesierr.png_masked 369193
## 23 23 MeltGrid.dat_masked 369193
## 24 24 mgsi.bmp_masked 369193
## 25 25 mgsierr.png_masked 369193
## 26 26 ssi.bmp_masked 369193
## 27 27 ssierr.png_masked 369193
## 28 28 subregions_Nittler_Vflip.bmp_masked 369193# ==================== SAUVEGARDE ====================
saveRDS(result_array_full, file = result_file)
cat("\n💾 Sauvegardé :", result_file, "\n")##
## 💾 Sauvegardé : /Users/alexandremichaux/Documents/UCA/Cours/Analyse des données/TP/TPs/Projet final/result_array_with_uncertainty.rdsresult_array_full[result_array_full == 0 | abs(result_array_full) < 1e-10] <- NA_real_
#nombre de NA dans chaque couche
na_counts <- sapply(1:dim(result_array_full)[3], function(i) sum(is.na(result_array_full[,,i])))
layer_info <- data.frame(Index = seq_along(attr(result_array_full, "layer_names")), Layer = attr(result_array_full, "layer_names"), NA_Count = na_counts)
print(layer_info) # Afficher les informations des couches avec le nombre de## Index Layer NA_Count
## 1 1 alsi.bmp 2237
## 2 2 alsierr.png 2236
## 3 3 casi.bmp 193559
## 4 4 casierr.png 194050
## 5 5 DEM.RDS 0
## 6 6 DensityGrid.dat 0
## 7 7 fesi.bmp 366605
## 8 8 fesierr.png 367271
## 9 9 MeltGrid.dat 0
## 10 10 mgsi.bmp 2238
## 11 11 mgsierr.png 2236
## 12 12 ssi.bmp 196990
## 13 13 ssierr.png 197926
## 14 14 subregions_Nittler_Vflip.bmp 811124
## 15 15 alsi.bmp_masked 369193
## 16 16 alsierr.png_masked 369193
## 17 17 casi.bmp_masked 369193
## 18 18 casierr.png_masked 369462
## 19 19 DEM.RDS_masked 369193
## 20 20 DensityGrid.dat_masked 369193
## 21 21 fesi.bmp_masked 369193
## 22 22 fesierr.png_masked 369859
## 23 23 MeltGrid.dat_masked 369193
## 24 24 mgsi.bmp_masked 369193
## 25 25 mgsierr.png_masked 369193
## 26 26 ssi.bmp_masked 369193
## 27 27 ssierr.png_masked 369363
## 28 28 subregions_Nittler_Vflip.bmp_masked 977473La fonction ci-dessous permet d’extraire une matrice 2D (une couche) du tableau 3D en fonction de l’index de la couche.
| Index | Layer |
|---|---|
| 1 | alsi.bmp |
| 2 | casi.bmp |
| 3 | CrustalThickness_ModelV1.csv |
| 4 | DEM.tif |
| 5 | DensityGrid.dat |
| 6 | fesi.bmp |
| 7 | MeltGrid.dat |
| 8 | mgsi.bmp |
| 9 | ssi.bmp |
| 10 | subregions_Nittler_Vflip.bmp |
get_layer_as_matrix <- function(result_array, layer_index) {
if (layer_index < 1 || layer_index > dim(result_array)[3]) { #check si l'index est inf. à 1 (index minimum) ou si il est sup. au nombre de couche de result_array
stop("Layer index out of bounds.")
}
return(result_array[,,layer_index]) # Extrait et retourne la couche spécifiée sous forme de matrice
}
for (i in 1:28){
layer_matrix <- get_layer_as_matrix(result_array_full, i)
image(
z = t(apply(layer_matrix, 2, rev)), # rotation correcte pour l'affichage
col = terrain.colors(100),
main = paste("Couche", i),
axes = FALSE
)}
Calcul de résidu en corrélation avec les expéreiences de fusion partielle….
# ==================== CALCUL DES CARTES DE PRESSION ET DE FUSION ====================
library(plotly)
# --- 1. Chargement des données expérimentales Mer8 et Mer15 ---
mer8_path <- "/Users/alexandremichaux/Documents/UCA/Cours/Analyse des données/TP/TPs/Projet final/data/data_Mer8.csv"
mer15_path <- "/Users/alexandremichaux/Documents/UCA/Cours/Analyse des données/TP/TPs/Projet final/data/data_Mer15.csv"
data_Mer8 <- read.csv(mer8_path, sep = ",", header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)
data_Mer15 <- read.csv(mer15_path, sep = ",", header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)
# Combiner les deux jeux expérimentaux
exp_data <- rbind(data_Mer8, data_Mer15)
exp_data <- exp_data[complete.cases(exp_data[, c("Mg/Si","Al/Si","Ca/Si","Fe/Si","S/Si")]), ]
# --- 2. Extraction des couches masquées depuis result_array_full ---
maps <- list(
MgSi = get_layer_as_matrix(result_array_full, 24),
MgSi_err = get_layer_as_matrix(result_array_full, 25),
AlSi = get_layer_as_matrix(result_array_full, 15),
AlSi_err = get_layer_as_matrix(result_array_full, 16),
CaSi = get_layer_as_matrix(result_array_full, 17),
CaSi_err = get_layer_as_matrix(result_array_full, 18),
FeSi = get_layer_as_matrix(result_array_full, 21),
FeSi_err = get_layer_as_matrix(result_array_full, 22),
SSi = get_layer_as_matrix(result_array_full, 26),
SSi_err = get_layer_as_matrix(result_array_full, 27)
)
nx <- nrow(maps$MgSi)
ny <- ncol(maps$MgSi)
# --- 3. Matrices vides pour Pression et F ---
pressure_map <- matrix(NA, nrow = nx, ncol = ny)
fusion_map <- matrix(NA, nrow = nx, ncol = ny)
# --- 4. Fonction de calcul du résidu pondéré ---
compute_residual <- function(M, sigma, E) {
valid <- !is.na(M) & !is.na(sigma)
if (sum(valid) < 3) return(NA)
sqrt(sum(((M[valid] - E[valid]) / sigma[valid])^2))
}
# --- 5. Boucle principale : minimisation du résidu pondéré ---
pb <- txtProgressBar(min = 0, max = nx, style = 3)## | | | 0%for (x in 1:nx) {
for (y in 1:ny) {
M <- c(maps$MgSi[x,y], maps$AlSi[x,y], maps$CaSi[x,y],
maps$FeSi[x,y], maps$SSi[x,y])
sigma <- c(maps$MgSi_err[x,y], maps$AlSi_err[x,y], maps$CaSi_err[x,y],
maps$FeSi_err[x,y], maps$SSi_err[x,y])
if (all(is.na(M))) next
residuals <- apply(exp_data[, c("Mg/Si","Al/Si","Ca/Si","Fe/Si","S/Si")], 1, function(E){
compute_residual(M, sigma, E)
})
best_index <- which.min(residuals)
if (is.na(best_index)) next
pressure_map[x,y] <- exp_data$Pression[best_index]
fusion_map[x,y] <- exp_data$F[best_index]
}
setTxtProgressBar(pb, x)
}## | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%close(pb)lst <- list(
pressure_map = pressure_map,
fusion_map = fusion_map
)
for (name in lst) {
image(
z = t(apply(name, 2, rev)), # rotation correcte pour l'affichage
col = terrain.colors(100),
main = "",
axes = T,
xlab = "Longitude",
ylab = "Latitude"
)}
# --- 7. Sauvegarde des résultats ---
saveRDS(pressure_map, file.path(base_path, "pressure_map_masked_weighted.rds"))
saveRDS(fusion_map, file.path(base_path, "fusion_map_masked_weighted.rds"))
cat("\n✅ Cartes de pression et de fusion pondérées sauvegardées.\n")##
## ✅ Cartes de pression et de fusion pondérées sauvegardées.Nittler, L.R., Frank, E.A., Weider, S.Z., Crapster-Pregont, E., Vorburger, A., Starr, R.D. & Solomon, S.C., 2020. global major-element maps of Mercury from four years of MESSENGER X-Ray Spectrometer observations. Icarus, 345, 113716. https://doi.org/10.1016/j.icarus.2020.113716.